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Preface 


The  objective  of  the  current  study  was  to  numerically 
solve  the  full  Nav i er -Stokes  equations  for  flow  over  an  iced 
NACA  8812  airfoil.  The  Beam-Warming  algorithm  was  used  to 
obtain  steady-state  solutions  for  angles  of  attack  ranging 
from  2-8  degrees.  Comparison  of  these  results  was  made  with 
other  numerical  results  obtained  from  algorithms  based  on 
approximate  forms  of  the  Navi er-Stokes  equations.  The 
present  results  were  compared  to  experimental  data  as  well. 

These  comparisons  revealed  that  the  Beam-Warming  scheme 
accurately  predicts  the  iced  airfoil  lift  and  drag  for 
angles  of  attack  below  stall.  The  post-stall  lift  is  over 
predicted,  however,  and  a  study  of  the  unsteady  nature  of 
the  flowfield  that  results  after  stall  will  be  needed  to 
resolve  this  problem.  Further  research  is  also  needed  to 
refine  the  computed  local  results,  such  as  the  reattachment 
location  of  the  separation  bubbles  aft  of  the  ice  shape. 
Continuation  of  this  effort  in  future  thesis  work  could  lead 
to  resolution  of  these  issues. 
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received  considerable  assistance  from  others.  I  would 
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research  effort.  From  start  to  finish,  Dr.  Halim  was  always 
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and  resolve  the  many  problems  that  I  faced  along  the  way.  I 
am  also  deeply  indebted  to  my  sponsor.  Dr.  Joe  Shang  of  the 
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the  computer  resources  used  in  this  project.  He  provided 
much  more  than  financial  support,  though.  Dr.  Shang  also 
took  time  from  his  busy  schedule  to  provide  very  valuable 
technical  advice,  and  for  his  help  in  this  respect  I  am 
quite  grateful.  !  also  wish  to  thank  Dr.  Miguel  Visbal  of 
the  Flight  Dynamics  Laboratory  for  his  assistance  with  the 
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Abstract 


The  accumulation  of  ice  on  an  aircraft  reduces  the  lift 
and  increases  the  drag  of  the  aircraft,  thus  posing  a 
serious  hazard  to  flight  safety.  The  primary  goal  of 
research  into  the  aircraft  icing  problem  is  to  develop 
methods  for  quantifying  the  detrimental  effects  on  aircraft 
performance.  This  task  can  be  handled  experimentally  or 
with  a  numerical  approach,  or  with  a  combination  of  both 
techniques.  Current  research  interest  focuses  on  developing 
computer  codes  which  can  be  used  to  predict  iced  airfoil 
performance.  The  experimental  data  gathered  over  the  years 
provides  a  valuable  database  for  verifying  the  codes  in 
deve I opmen  t . 

The  objective  of  the  present  study  is  to  evaluate  the 

performance  of  an  iced  NACA  8812  airfoil  numerically.  The 

full  Nav i er-Stokes  equations  are  solved  using  the 

Beam-Warming  implicit  factored  scheme.  Steady-state 

solutions  are  obtained  for  flow  at  a  free  stream  Mach  number 

0 

of  8.12  and  a  Reynolds  number  based  on  chord  of  1.41X18  , 
for  angles  of  attack  of  2,  4,  6,  and  8  degrees. 

The  lift  and  drag  curves  obtained  from  the  numerical 
solutions  were  compared  to  experimental  data  and  other 
numerical  results.  It  was  found  that  the  Beam-Warming 
algorithm  provides  a  good  estimate  of  the  lift  and  drag  at 
angles  of  attack  below  stall.  Computed  lift  coefficients 


I 


were  within  11.5%  of  the  experimental  data,  while  the  drag 
coefficients  differed  by  as  much  as  24%.  These  results  were 
in  excellent  agreement  with  the  other  numerical  solutions. 

At  an  angle  of  attack  above  stall,  however,  the  code  did  not 
predict  the  expected  decrease  in  lift,  and  the  calculated 
drag  coefficient  was  much  lower  than  the  experimental  data. 

The  comparison  of  two  local  flowfield  char ac ter i st i cs 
with  the  experimental  data  was  less  encouraging.  A  computed 
velocity  profile  was  compared  to  the  experimental  profile 
for  a  station  in  the  separated  region  on  the  upper  surface. 
This  comparison  showed  that  the  computed  separation  bubble 
is  approximately  one  half  the  thickness  of  the  bubble 
measured  experimentally.  Flow  reattachment  location  is 
another  measure  of  the  predict ive  accuracy  of  the  numerical 
scheme.  The  reattachment  locations  from  the  numerical 
solutions  were  38-48%  less  than  the  experimental  values. 

The  conclusion  to  be  drawn  from  these  comparisons  is 
that  accurate  results  can  be  obtained  for  the  global 
performance  parameters  using  the  code  in  its  present  state, 
but  more  work  is  certainly  needed  to  refine  the  local 
results.  To  achieve  this  goal,  further  study  is  needed  in 
three  areas:  1)  effect  of  the  grid  on  solution  accuracy,  2) 
effect  of  the  transition  location  on  the  results,  and  3) 
effect  of  the  turbulence  model  on  the  solution. 
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NUMERICAL  SIMULATION  OF  FLOW  OVER  ICED  AIRFOILS 

I  .  I n t roduc t i on 

Ice  accumulation  on  aircraft  surfaces  results  in 
substantial  performance  penalties.  Formation  of  ice  on  the 
leading  edge  of  an  airfoil  can  lead  to  separated  flow 
regions  on  the  upper  and  lower  surfaces.  Flow  separation 
results  in  severe  loss  of  lift  and  a  dramatic  increase  in 
the  drag  of  the  airfoil.  For  example,  in  experimental  work 
at  Ohio  State  University,  Bragg  has  tested  a  NACA  0012 
airfoil  with  a  simulated  ice  shape  which  resulted  in  a 
reduction  of  maximum  lift  of  over  50%  and  a  drag  increase  of 
300%  (1:1).  The  implications  of  such  severe  performance 
penalties  on  flight  safety  are  obvious. 

An  aircraft  can  encounter  icing  when  flying  through 
clouds.  Two  types  of  ice  can  form  on  the  aircraft  surfaces 
depending  upon  the  environmental  conditions  (2:2-3).  When 
the  air  temperature  is  well  below  freezing  and  the  liquid 
water  content  of  the  air  is  low,  the.  water  droplets  freeze 
almost  immediately  upon  impact  forming  spherical  grains  of 
ice.  This  type  of  ice  formation  is  termed  rime  ice.  If  the 

air  temperature  is  near  freezing  and  the  liquid  water 
content  of  the  air  is  high,  the  droplets  impinging  on  the 
surface  do  not  freeze  immediately.  Instead,  droplets  run 
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back  on  the  airfoil  surface  before  freezing  which  leads  to 
the  formation  of  horn-shaped  structures  on  the  upper  and 
lower  airfoil  surfaces.  This  type  of  ice  formation  is 
called  glaze  ice.  The  glaze  ice  horns  cause  the  flow  to 
separate.  Therefore,  glaze  ice  is  the  most  detrimental  to 
aircraft  performance.  Typical  rime  and  glaze  ice  formations 
are  shown  in  Figure  1. 

The  problem  of  aircraft  icing  has  been  of  interest  for 
many  years.  In  1932,  Jacobs  tested  a  NACA  QQ 1 2  airfoil  with 
various  leading  edge  protrusions  (3).  Although  none  of 
these  experiments  were  conducted  with  actual  ice  accretions, 
the  detrimental  effect  of  distorting  the  leading  edge  shape 
was  established  in  this  work.  Jacobs’  results  were  reported 
in  terms  of  changes  in  the  airfoil's  lift,  drag,  and  moment. 

In  the  1956s  Gray  continued  NACA’s  study  with  tests  of 
the  NACA  65-212  and  NACA  65A664  airfoils.  These  experiments 
were  conducted  in  NACA’s  Icing  Research  Tunnel  and  data  were 
obtained  for  these  airfoils  with  actual  ice  accretions.  The 
adverse  effect  of  the  ice  was  again  reported  in  terms  of  the 
global  airfoil  performance  par ameter s-- I i f t ,  drag,  and 
moment  coefficients  (4,5).  In  1964,  Gray  used  the 
experimental  data  available  at  the  time  to  correlate  changes 
in  aerodynamic  performance  with  ice  accumulation.  He  also 
attempted  to  generalize  the  correlation  to  other  airfoils 
using  leading  edge  radius  as  a  parameter  (6:1).  While  this 
correlation  has  proven  useful  in  a  qualitative  sense,  Cebec i 


points  out  that  the  drag  rise  predicted  with  Gray's 
correlation  is  generally  too  high  (7:11). 

The  need  to  more  accurately  predict  the  effect  of  ice 
accumulation  on  airfoils  has  prompted  further  research  in 
this  area  in  the  1989s.  The  ultimate  goal  of  the  current 
NASA  program,  for  instance,  is  to  develop  computer  codes 
which  can  be  used  in  the  airfoil  design  process  to  predict 
the  sensitivity  of  a  proposed  design  to  ice  accumulation. 
Before  any  code  can  be  used  as  a  reliable  design  tool, 
however,  it  must  be  thoroughly  tested.  One  method  to  test  a 
code  under  development  is  to  apply  the  code  to  an  iced 
airfoil  for  which  experimental  data  are  available,  and 
compare  airfoil  performance  based  on  the  computed  solution 
to  that  obtained  from  the  experimental  data.  Obviously, 
accurate,  detailed  experimental  data  for  an  existing  iced 
airfoil  are  required  to  fully  validate  a  given  code. 

In  recent  years,  several  investigators  have  worked  to 
provide  the  detailed  experimental  data  needed  to  verify  the 
computer  codes  in  development.  Bragg  and  Coirier  have 
performed  a  series  of  experiments  using  a  NACA  8912  airfoil 
with  a  simulated  glaze  ice  shape  attached  to  the  leading 
edge.  In  these  experiments,  the  lift,  drag,  and  moment  of 
the  airfoil  was  determined  for  the  airfoil  with  and  without 
ice.  In  addition,  surface  pressure  taps  were  located  in  the 
regions  behind  the  ice  horns  to  allow  measurement  of  the 
pressure  distribution  in  the  separated  flow  regions.  Also, 
a  split  film  probe  was  used  to  map  out  the  velocity  profiles 
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at  several  chordwise  stations  in  the  separation  bubble 
behind  the  ice  horns.  These  velocity  profiles  clearly 
indicate  the  regions  of  reversed  flow  in  the  bubble  and  the 
bubble  reattachment  point  (1,8).  In  an  extension  of  this 
work,  Bragg  and  Spring  have  investigated  the  effect  of 
adding  surface  roughness  to  the  glaze  ice  shape  (9,18).  in 
addition  to  these  quantitative  studies,  Khodadoust's  flow 
visualization  study  of  this  same  airfoil  provides  some 
qualitative  information  on  the  behavior  of  the  flow  in  the 
leading  edge  separation  bubble  (11). 

In  addition  to  the  experimental  work  on  iced  airfoils 
in  recent  years,  there  has  been  an  effort  to  develop 
computational  methods  for  determining  the  flow  field  over 
iced  airfoils.  Two  researchers  ,  Potapczuk  at  NASA’s  Lewis 
Research  Center  and  Cebeci  at  the  California  State 
University,  have  published  very  promising  results  for  the 
flow  over  the  NACA  8812  airfoil  tested  experimentally  by 
Bragg  and  Coirier.  Both  investigators  have  produced 
computational  results  which  compare  very  favorably  to  the 
lift  and  drag  data  of  Bragg  and  Coirier,  even  though 
different  numerical  schemes  are  used  by  each. 

Cebeci  has  developed  a  code  for  iced  airfoils  based  on 
an  interactive  boundary  layer  (IBL)  method  (12).  The  code 
allows  for  interaction  of  the  viscous  and  inviscid  flow 
regimes.  With  this  feature,  the  boundary  layer  equations 
are  well  behaved  even  beyond  the  point  of  separation  on  the 
airfoil  surface.  Therefore,  the  code  can  be  marched  past 


the  point  of  separation  and  a  solution  can  be  obtained  even 
in  regions  of  reverse  flow.  Turbulent  eddy  viscosity  is 
computed  using  Cebeci  and  Smith's  algebraic  turbulence 
model.  Details  of  the  IBL  method  are  in  Reference  12.  The 
Cebeci-Smith  turbulence  model  is  described  in  Reference  13. 

Potapczuk  developed  a  computational  method  for  the  flow 
over  iced  airfoils  using  the  thin  layer  Nav i er -St okes 
equations  (14,15).  The  code  being  used  is  ARC2D,  a 
Nav i er-Stokes  code  developed  by  Steger  and  Pulliam  at  NASA 
Ames.  The  algebraic  turbulence  model  of  Baldwin  and  Lomax 
is  used  in  this  code  to  calculate  the  turbulent  eddy 
viscosity.  Details  of  the  ARC2D  code  can  be  found  in 
References  16  and  17.  The  Ba I dwi n-Lomax  turbulence  model  is 
described  fully  in  Reference  18. 

In  the  present  work,  two  additional  numerical  methods 
are  investigated  for  calculating  the  flow  field  over  the 
NACA  8812  airfoil  with  glaze  ice.  The  goal  of  the  study  is 
to  compare  the  present  numerical  results  to  other  numerical 
results  and  experimental  data,  and  evaluate  the  advantages 
or  disadvantages  of  the  algorithms  studied  in  this  thesis. 

To  match  the  experimental  conditions  of  Bragg  and  Coirier, 
the  numerical  solutions  are  obtained  for  a  Mach  number  of 
8.12  and  a  Reynolds  number  of  1.41X18  .  Numerical  solutions 
are  obtained  for  ang I es-o f -a t t ack  ranging  from  2-8 
degrees . 

One  of  the  numerical  methods  investigated  is  a  global 
space  marching  technique  developed  by  Halim  (19).  This  code 


is  based  on  an  approximate  form  of  the  Nav i er-Stokes 
equations  obtained  by  dropping  the  streamwise  diffusion 
terms.  The  equations  are  written  in  terms  of  the  stream 
function  and  vorticity  rather  than  the  primitive  variables. 
The  advantage  of  this  formulation  is  that  the  number  of 
unknowns  is  reduced  by  one. 

A  second  numerical  method  applied  to  the  iced  airfoil 
is  the  implicit  approximate  factored  scheme  due  to  Beam  and 
Warming  (28).  This  code  is  a  time  marching  implementation 
of  the  full  Nav i er-Stokes  equations.  No  simplifying 
assumptions  have  been  made.  This  code  is  written  in  terms 
of  the  primitive  variables. 

There  are  two  common  features  for  these  codes.  First 
of  all,  each  code  is  written  in  terms  of  generalized 
coordinates  to  allow  computation  over  arbitrary  geometries. 
This  is  accomplished  by  transforming  an  arbirtray  physical 
space  to  a  computational  space  with  uniform  spacing. 
Secondly,  each  code  is  written  in  delta  formulation.  This 
means  that  the  solution  obtained  is  in  terms  of  the 
difference  between  a  variable’s  value  at  two  time  levels. 

In  the  limit,  this  difference  approaches  zero  as  the 
solution  approaches  steady-state. 


* 


I  I .  Ana  lysis 

* 

In  this  chapter  the  governing  equations  for  flow  over 
an  airfoil  are  described.  Anderson  and  Tannehill  (21)  is 
*  the  primary  reference  for  the  equations  presented  in  this 

chapter.  The  governing  equations  are  put  into 
non-dimensional  form  and  transformed  from  Cartesian 
1  coordinates  to  generalized  coordinates.  The  method  of 

generating  the  grid  for  solving  the  equations  in  the 
transformed  domain  is  briefly  discussed.  Finally,  the 
’  appropriate  boundary  conditions  for  flow  over  an  airfoil  are 

presented . 

1  Govern i ng  Equa  t i on s 

The  flow  field  over  an  iced  airfoil  is  described  by  the 
Navier-Stokes  equations.  The  vector  form  of  the 

i 

compressible  Navier-Stokes  equations  in  Cartesian 
coordinates  is  given  by 


dU  3E  ^  dF  0G  _ 

3F  +  3*  +  5y  +  57  - 

These  equations  can  be  simplified  for 
t  o 


Q 


two-d i mens i ona  I 


(  1 ) 

f  I  ow 


au  be  A  aF 
3T  +  3x  +  ay 


s 


(2) 


8 


where  the  vector  of  conserved  variables,  U,  is  given  by 
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and  the  flux  vectors,  E  and  F  are 
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For  turbulent  flow  the  shear  stress  components  are  expressed 
as 
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where 

/i  =  molecular  dynamic  viscosity 
e  =  turbulent  eddy  viscosity 

Also,  assuming  Stokes’  hypothesis  to  be  valid  for  this  flow 
Xj  is  def i ned  as 


XT  =  -  2/3  {  fi  +  e  ) 


The  heat  flux  components  for  turbulent  flow  are 


(  k  *  kT  ) 
(  k  +  kT  ) 


where 


thermal  conductivity  of  air 
turbulent  thermal  conductivity 


Introducing  the  Prandtl  number,  Pr ,  and  turbulent  Prandtl 
number,  Pfj,  defined  as 


Pr  = 


PrT  = 


the  heat  flux  terms  can  be  rewritten  as 
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The  energy  term,  E{,  in  the  E  and  F  vectors  above  denotes 
the  total  energy  per  unit  volume,  given  by 


E*  =  P 
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where 

u,  v  =  velocity  components 

e  =  internal  energy  per  unit  mass 

The  equation  of  state  for  a  perfect  gas  is  used  to  close  the 
system  of  equations.  This  equation  is 


P  =  pRT  (16) 

where  R  is  the  gas  constant  for  air.  Using  the  relation  for 
a  calorically  perfect  gas 
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where 
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cp  =  specific  heat  at  constant  pressure 
cy  =  specific  heat  at  constant  volume 

the  internal  energy  can  be  expressed  as 


e 


(28) 


Substituting  for  RT  in  the  equation  of  state  yields 


P  =  pe  (  r  -  1  ) 


(21) 


Using  Eq  (21)  the  pressure  can  be  expressed  in  terms  of  the 
total  energy  as 
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In  addition  to  the  equations  developed  above, 
expressions  are  needed  to  evaluate  the  molecular  and  eddy 
viscosities.  The  molecular  viscosity  is  determined  using 
Sutherland's  formula  (22:328),  given  by 
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where 

Tq  r  absolute  reference  temperature 
p  =  molecular  viscosity  at  reference  temperature 
S1  s  Suthe'land’s  constant  {  118  K  ) 
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The  turbulent  eddy  viscosity  is  evaluated  using  the 
algebraic  eddy  viscosity  model  of  Baldwin  and  Lomax  (18),  as 
modified  by  Visbal  (23). 


Non-dimens i ona I  Form  of  the  Equa t i ons 

It  is  customary  to  cast  the  governing  equations  into 
non-dimensional  form  so  they  are  expressed  in  terms  of 
parameters  such  as  Reynolds  number,  Mach  number  and  Prandtl 
number.  One  method  for  non-d imens i ona I i 2 i ng  the  equations, 
as  suggested  by  Anderson  and  Tannehill  (21),  is  to  define 
dimensionless  variables  of  the  form 


=  x/c 

y*  =  y/c 

* 

t  -  - 

t 

c7V 

CO 

=  u/V 

OD 

v*  =  v/v 

OD 

* 

p  = 

=  ^QD 

£*  =  e/*l„ 

* 

T  = 

T/T» 

* 

p  = 

P  .  e*  = 

e 

2  e 
p\r 

K  00 

Vs 

00 

(24) 


where 


c  =  airfoil  chord  length 
V  =  free  stream  velocity 

OD 


The  asterisk  denotes  the  dimensionless  variables  and  the 
subscript  00  indicates  free  stream  conditions.  Applying 
these  expressions  to  the  Nav i er-Stokes  equations  yields 
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where  the  dimensionless  vectors  are  expressed  as 
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In  these  vectors  the  dimensionless  shear  stress  terms  are 


given  by 
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The  corresponding  heat  flux  expressions  are 
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The  dimensionless  form  of  the  total  energy  is  expressed  as 
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and  the  corresponding  dimensionless  equations  of  state  are 
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Note  that  the  dimensionless  vectors  have  exactly  the  same 
form  as  the  dimensional  ones.  For  this  reason,  the  * 
notation  is  generally  dropped  once  the  dimensionless  form  of 
the  equations  has  been  derived.  This  convention  will  be 
followed  in  this  thesis.  In  all  subsequent  references  to 
the  governing  equations  the  dimensionless  form  is  implied. 
Details  of  the  non-d imens i ona I i zat i on  for  representative 
terms  are  shown  in  Appendix  A. 


Equations  in  Generalized  Coordinates 


The  equations  as  developed  so  far  are  suitable  for 
application  to  problems  in  the  Cartesian  coordinate  system. 
For  problems  with  flow  over  irregular  geometries,  however, 
the  equations  can  be  more  conveniently  solved  if  expressed 
in  terms  of  generalized  coordinates.  The  general 
transformation  given  by 


5  =  5(x,y) 

r?  =  r?(  x  ,  y  ) 


can  be  applied  to  transform  the  governing  equations  from  the 
physical  domain  { x ,  y )  to  a  computational  domain  (£,n)-  The 
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transformation  is  shown  in  Figure  2.  As  can  be  seen,  the 
advantage  of  the  transf ormat i on  is  that  a  uniformly  spaced 
computational  grid  results. 

As  outlined  in  Anderson  and  Tannehill  (21),  the 
transf ormat i on  from  physical  coordinates  to  computational 
coordinates  is  performed  by  applying  the  chain  rule  of 
partial  differentiation  as  follows: 
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where  J  is  the  Jacobian  of  the  transformation  defined  as 
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Further  details  of  the  coordinate  transformation  are  in 


Appendix  B. 

Now  Eq  (41)  can  be  rewritten  as 


where 
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One  further  rearrangement  of  the  governing  equations  is 
needed  to  aid  in  the  implementation  of  the  implicit 
algorithm.  It  is  desirable  to  separate  the  inviscid  and 
viscous  terms  in  the  E  and  F  vectors,  and  it  is  also 
necessary  to  separate  the  cross-derivative  terms.  This  is 
accomplished  by  rearranging  terms  into  six  new  vectors. 

The  inviscid  and  pressure  terms  in  E  in  Eq  (44)  form  a 
new  vector  E^.  Similarly,  the  inviscid  and  pressure  terms 
in  F  form  a  new  vector  Eg.  These  vectors  can  be  expressed 
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where  the  cont r avar i ant  velocities,  U  and  V,  are  defined  as 
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The  viscous  and  heat  flux  terms  which  remain  in  Eq  (44) 
are  separated  into  four  vectors  to  isolate  the 
cr oss-der i va t i ves .  These  terms  in  E  can  be  expressed  in 
general  as 


V  =  V,(G,U5)  *  V^G.uJ  (52, 

and  the  terms  in  F  can  be  written  in  general  as 

W  =  W^U.uJ  +  W^U.uJ  (53) 


With  these  definitions,  Eq  (44)  can  be  rewritten  as 


The  specific  components  of  V^,  W1  ,  and  involve  the 

expressions  for  the  shear  stress  and  heat  flux  terms  in 
generalized  coordinates.  These  expressions  are  included  in 
Appendix  B. 

The  Navier-Stokes  equations  are  now  in  a  form  suitable 
for  application  to  the  Beam-Warming  implicit  algorithm. 
Details  of  this  numerical  scheme  are  discussed  in  the  next 
chapter. 

Gr i d  Gener at i on 

In  the  numerical  study  of  any  fluid  dynamics  problem, 
the  continuous  solution  is  approximated  by  solving  the 
governing  equations  at  discrete  points  in  the  flow  field. 
The  distribution  of  these  points  in  space  forms  a  grid  or 
mesh  system.  Accuracy  of  the  solution  for  a  given  problem 
is  highly  dependent  on  the  number  and  placement  of  the  grid 
points. 


For  simple  geometries,  the  grid  can  be  defined  by 
simply  specifying  the  spacing  and  number  of  points  in  the 
directions  parallel  and  normal  to  the  body.  The  governing 
equations  are  then  solved  directly  in  the  physical  domain. 
For  more  complicated  shapes,  such  as  two-dimensional 
airfoils,  fitting  the  grid  to  the  body  results  in  rather 
arbitrary  spacing  in  the  physical  domain.  Specifying  the 
boundary  conditions  on  the  surface  in  physical  coordinates 
can  also  be  difficult  for  complex  shapes.  Therefore,  the 
physical  domain  is  generally  transformed  into  a 
computational  domain.  The  transformation  is  such  that 
points  in  the  physical  domain  are  mapped  one-for-one  into  a 
unit  square  computational  space,  as  shown  in  Figure  2.  Note 
that  in  the  computational  domain  the  spacing  is  uniform 
throughout.  Also,  the  airfoil  surface  maps  to  the  q=B 
line  in  the  computational  domain  which  greatly  simplifies 
specification  of  the  surface  boundary  conditions. 

Several  techniques  are  available  for  performing  this 
transformation,  including  algebraic  methods  and  methods 
requiring  the  solution  of  elliptic  or  hyperbolic  partial 
differential  equations.  Each  has  advantages  and 
limitations.  In  this  study,  airfoil  grids  were  generated 
using  both  partial  differential  equation  methods. 

initial  work  was  done  using  an  elliptic  grid  generator 
developed  by  David  Amdahl  (24)  of  the  Flight  Dynamics 
Laboratory  at  Wr i ght-Pat terson  AFB ,  Ohio.  The  code  worked 
very  well,  but  the  number  of  iterations  required  to  solve 
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the  elliptic  equations  translated  into  long  computer  running 
times  for  each  grid.  Also,  the  version  of  the  code  used  in 
this  study  requires  the  user  to  specify  the  outer  boundary 
points  in  addition  to  the  inner  boundary  or  airfoil  surface 
points.  Given  these  limitations,  an  alternative  grid 
generation  method  was  chosen. 

Subsequent  grids  were  generated  using  a  hyperbolic  grid 
generator  developed  by  Steger  and  Chaussee,  and  coded  by 
Kinsey  and  Barth  (25).  The  solution  of  the  hyperbolic 
system  of  equations  used  in  this  code  was  obtained  much 
faster  than  with  the  elliptic  grid  generator.  Also,  since 
the  hyperbolic  equations  are  marched  in  space,  only  an  inner 
boundary  must  be  specified  by  the  user.  The  program  allows 
for  selection  and  modification  of  the  desired  outer  boundary 
shape.  For  example,  C-shaped  or  O-shaped  grids  can  be 
generated  about  an  airfoil  simply  by  adjusting  program 
default  parameters.  In  the  present  study,  C-grids  were 
generated.  A  grid  with  299  points  in  the  ^-direction  and  65 
points  in  the  redirection  is  shown  in  Figure  3. 

The  hyperbolic  grid  generator  offers  complete  control 
over  grid  point  clustering  parallel  and  normal  to  the 
surface.  For  the  iced  airfoil,  it  is  necessary  to  cluster 
more  points  around  the  leading  edge  region  to  capture  the 
rapid  changes  in  the  flow  over  the  ice  shape.  The  grid 
generator  allows  the  user  to  specify  the  minimum  spacing 
along  the  surface  at  any  point,  either  on  the  body  or  in  the 
wake.  Once  these  points  are  all  defined,  the  user  can 
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specify  the  number  of  points  needed  between  any  two  minimum 
spacing  locations.  Exponential  stretching  functions  are 
used  to  space  these  grid  points  on  the  surface  or  in  the 
wake  about  the  minimum  spacing  points  previously  defined. 

Similarly,  the  code  allows  the  user  to  specify  the 
minimum  spacing  required  normal  to  the  airfoil  surface. 

This  provides  the  control  needed  to  cluster  more  points  in 
the  boundary  layer  where  rapid  flow  changes  are  occurring. 
Again,  exponential  stretching  is  applied  to  space  the 
remaining  points  normal  to  the  surface. 

Boundary  Cond i t i ons 

For  solution  of  the  compressible  Nav i er-Stokes 
equations,  boundary  conditions  must  be  specified  at  the  grid 
outer  boundary  and  on  the  airfoil  surface.  It  is  assumed 
that  the  outer  boundary  is  far  enough  away  from  the  airfoil 
that  free  stream  conditions  can  be  applied.  Therefore,  in 
the  computational  domain  at  rp n_._v  .  the  boundary 

cond i t i ons  are 
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where 


00 


=  magnitude  of  the  free  stream  velocity 


a  =  airfoil  angle  of  attack 


At  the  downstream  boundaries,  §=8  and  €=£_av 
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specified  boundary  conditions  are 
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For  the  viscous  flow  on  the  airfoil  surface,  the  no-slip 
boundary  condition  is  enforced.  An  adiabatic  wall  is 
assumed,  and  the  pressure  at  the  surface  is  assumed  to  be 
equal  to  the  pressure  one  grid  point  above  the  surface. 
Thus,  for  rp8  ,  the  boundary  conditions  are 
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(57) 


I  In  the  wake,  continuity  must  be  maintained  across  the  wake 

cut.  This  is  accomplished  by  averaging  values  of  p,  u,  v, 
and  p  across  the  wake  cut  and  using  these  average  values  as 
>  the  boundary  condition  on  the  wake  line. 


i 


111.  Numerical  Procedures 


The  numerical  methods  used  to  solve  the  governing 
equations  are  described  in  this  chapter.  Two  numerical 
schemes  have  been  used  to  study  the  flow  over  a  NACA  8812 
a i r f o i I -- the  Beam-Warming  finite-difference  scheme  (28,26) 
and  a  space  marching  algorithm  developed  by  Halim  (19).  The 
Beam-Warming  scheme  is  an  implementation  of  the  full 
compressible  Nav i er-Stokes  equations.  The  actual  code  used 
was  written  by  Visbal  (23).  The  Halim  formulation  is  an 
implementation  of  the  incompressible  Nav i er-Stokes  equations 
utilizing  the  thin  shear  layer  approximation.  It  was 
originally  written  and  applied  to  symmetric  airfoils  at  zero 
angle  of  attack.  By  taking  advantage  of  the  symmetry 
condition,  the  solution  only  needed  to  be  computed  over  half 
the  domain.  The  author  and  a  fellow  AFIT  student,  Capt 
Faran  Hafeez  (27),  worked  to  modify  the  code  to  handle 
computations  in  the  full  domain. 

Beam-Warm i ng  Scheme 

The  Beam-Warming  scheme  is  a  time  marching  implicit 
algorithm  for  the  solution  of  the  compressible  Nav i er -S t okes 
equations.  The  general  difference  formula  for  the  scheme  is 
given  by 
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(58) 


where 

AUn  =  Un+1  -  Un  (59) 

The  temporal  accuracy  of  the  scheme  is  determined  by  the 
choice  of  the  parameters  fl  ,  and  0 2>  The  Euler  implicit 
scheme  is  obtained  for  0^  =  1  and  e2  =  8  .  This  scheme 
is  first  order  accurate  in  time.  A  variety  of  second  order 
accurate  schemes  can  be  obtained  by  specifying 
01  =  02  +  1/2  .  The  specific  code  used  in  this  thesis  uses 

first  order  accurate  Euler  implicit  time  differencing.  This 
accuracy  is  sufficient  in  the  present  work  since  only  the 
steady-state  solutions  are  of  interest. 

Recall  that  the  Nav i er-Stokes  equations  are  given  by 


(54) 
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If  this  equation  is  written  at  time  levels  n+1  and  n 
and  the  equation  at  level  n  is  subtracted  from  the 
equation  at  level  n+1  ,the  result  is 
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This  is  the  so-called  delta  form  of  the  equations,  where 
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Now  Eq  (68)  can  be  solved  for  AUt  and  the  result  substituted 
into  Eq  ( 58 )  to  yield 
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For  the  Beam-Warming  algorithm,  the  spatial  derivatives  in 
Eq  (62)  cannot  be  approximated  by  finite  differences  yet 
because  the  equation  is  nonlinear.  This  arises  from  the 
fact  that  the  flux  vector  increments,  AE^  AEg,  AV^  AVg, 
AW1 ,  and  AWg,  are  nonlinear  functions  of  U.  Each  of  the 
flux  vectors  must  be  linearized  before  differencing  the 
equa  t i on . 

Linear i zat i on  of  the  F I ux  Vectors.  Each  of  the  flux 
vectors  can  be  linearized  using  a  truncated  Taylor  series 
expansion  in  U.  For  example,  AE*J  can  be  linearized  from 


En+1  =  En  + 
C1  C1 
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which  can  be  expressed  as 


n  9 
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where  f  A  1  is  the  Jacobian  matrix  given  by 
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Similarly,  AEg  is  linearized  as 
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where 


s  the  Jacobian  matrix  given  by 


M' 


The  viscous  delta  terms  are  linearized  using  the  method 
suggested  by  Steger  (28).  In  this  method,  the  coefficients 
of  viscosity,  jj  and  e,  and  the  thermal  conductivities,  k  and 

A 

kT ,  are  assumed  to  be  locally  independent  of  U.  In 

addition,  the  cross  derivative  terms,  AV»  |U,U  ]  and 

2gl  nj 

Aw'J  ^U.U^j  ,  are  neglected.  With  these  approximations,  the 
viscous  terms  are  linearized  as 
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where  f  R  J  and  £  S  J  are  the  Jacobian  matrices 


(  7® 
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The  four  Jacobian  matrices  resulting  from  the  linearization 
of  the  governing  equations  are  presented  in  Appendix  C. 

Substituting  the  linearized  flux  vectors  into  Eq  (62) 
yields 
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where 

£  I  j  s  identity  matrix 
and  notat i on  such  as 
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imp  I i es 


Approximate  Factor i zat i on .  Beam  and  WarminQ  (28)  have 
shown  that  spatial  factorization  of  Eq  (72)  results  in  a  set 
of  equations  that  can  be  solved  using  an  alternating 
direction  numerical  scheme.  Specifically,  the  left  hand 
side  ( LHS )  OF  Eq  (72)  can  be  factored  into 
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Since  the  term  remaining  after  factorization  is  O(At)  ,  and 

the  leading  truncation  error  term  in  Eq  (72)  is  at  best 
2 

O(At)  ,  the  factorization  does  not  affect  the  temporal 
accuracy  of  the  scheme. 

In  practice,  Eq  (72)  is  implemented  in  three  steps  as 
foil ows : 

Step  1 . 
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=  RHS  Eq  (72) 


Step  3 . 

Un+1  -  un  +  Ayn  (7 

where 

... 

AU  =  intermediate  solution  vector 

RHS  £  Eq  (72)  j  =  terms  on  right  hand  side  of  Eq  (72) 

The  equations  above  are  written  in  terms  of  a  general 
time  differencing  algorithm.  The  code  used  in  this  thesis 
incorporated  first-order  Euler  implicit  time  differencing, 
obtained  by  setting  01  =  1  and  =  8  .  With  these 

substitutions  for  01  and  02 ,  Eqs  (74)  and  (75)  simplify  to 


(78] 


The  spatial  derivatives  in  Eqs  (77)  and  (78)  are  written  as 
second-order  central  differences.  Visbal  (23)  shows  that 
this  yields 
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are  finite  difference  operators. 


This  discretization  of  the  equations  results  in  a  4X4 
block  tridiagonal  system  of  linear  equations  to  be  solved. 
The  solution  is  obtained  by  first  sweeping  the  domain  in  the 
(►-direction,  solving  Eq  (79)  at  each  17-line.  This  sweep  is 
followed  by  sweeping  the  domain  in  the  rj-direction,  solving 
Eq  (86)  along  each  g-line.  The  result  of  both  sweeps  is  a 
new  estimate  for  the  unknowns,  AUn. 

Added  Dissipation  Terms .  In  the  application  of  their 
algorithm,  Beam  and  Warming  have  found  it  necessary  to  add 
artificial  dissipation  terms  to  damp  short  wave  length 
oscillations  (26).  In  Visbal's  implementation  of  the 
Beam-Warming  scheme,  both  explicit  fourth-order  damping 
terms  and  implicit  second-order  damping  terms  have  been 
added  for  numerical  stability.  Explicit  damping  is  included 
by  adding  the  term 

K  *  *i)  G".i  1821 

to  the  RHS  of  Eq  (79),  while  implicit  damping  is 
incorporated  by  including  the  terms 

-  1  ]  |S3> 

and 


(84) 


on  the  LHS  of  Eqs  (79)  and  (89).  In  these  expressions,  Ug 
and  <■>£  are  the  explicit  and  implicit  damping  factors, 
respectively.  Anderson  and  Tannehi I  I  (21)  suggest  <jg  be 
less  than  9.9625  for  stability,  and  Visbal  (23)  and  Pulliam 
(17)  recommend  that 

(j^  2  2b)g  (85) 

Note  that  the  addition  of  the  second-order  damping 
terms  does  not  affect  the  steady-state  solution,  since  AUn 
approaches  zero  for  this  condition. 

Hal im*  s  Formu I  a  t i on  o  f  t  he  Th i n-Layer  Equations 

Another  numerical  method  considered  in  the  present 
study  was  that  developed  by  Halim  (19)  for  steady, 
incompressible,  laminar  flow  over  arbitrary  airfoils.  This 
code  employs  the  thin-layer  form  of  the  Nav ier-Stokes .  The 
thin-layer  approximation  assumes  the  viscous  shear  stresses 
parallel  to  the  body  are  negligible  compared  to  the  shear 
stress  normal  to  the  body.  Therefore,  the  streamwise 
diffusion  terms  are  dropped  from  the  Nav i er -S t okes 
equa  t i ons . 

Halim’s  formulation  of  the  thin-layer  equations  is 
written  in  terms  of  the  stream  function  and  vorticity  rather 
than  the  primitive  variables.  The  advantage  of  this 
formulation  for  2-D  flow  is  that  the  number  of  unknowns  is 
reduced  by  one.  In  terms  of  stream  function  and  vorticity, 


the  governing  equations  can  be  written  in  generalized 
coordinates  as 
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where  the  metrics  and  Jacobian  are  given  by 
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The  vorticity,  <o,  is  defined  as 
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and  the  stream  function  is  defined  by  the  equations 
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Linearization  of  the  Vorticit 


i on .  Equation  (86) 


is  nonlinear  due  to  the  products  up.,  and  <ju  This 

I  n 

equation  can  be  linearized  by  writing 
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n+1  n.n  , _  _ . 

tp  =  tp  +  Atp  (93) 

and  substituting  these  expressions  into  Eq  (86).  Expanding 
after  the  substitution,  and  neglecting  higher  order  terms, 
resu Its  in 
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The  equation  is  now  linear  in  the  unknowns  Aw  and  A«p. 

Equation  (87)  is  already  linear,  but  it  can  be  put  in 
delta  form  by  applying  Eqs  (92)  and  (93),  which  gives 
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Equations  (94)  and  (95)  are  discretized 
central  differences  for  all  terms  except 
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tp  Aw_ 


second-order 
,  which  is 
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treated  as  an  upwind  difference  for  stability.  A  2X2  block 
tridiagonal  system  of  equations  results,  which  is  solved 
using  the  Thomas  algorithm  (21:99) 


Boundary  Cond i t i ons  i n  Terms  of  Stream  Func t i on  and 
Vor t i c i ty .  The  boundary  conditions  used  for  this  algorithm 
are  similar  to  those  applied  for  the  Beam-Warming  scheme. 
However,  the  appropriate  boundary  conditions  must  now  be 
specified  implicitly  in  terms  of  p  and  w. 

On  the  airfoil  surface,  the  velocity  components  are 
zero.  Specifying  this  condition  in  terms  of  the  stream 
function  requires  applying  the  chain  rule  to  the  stream 
function  definition,  which  yields 


u  =  Vy  =  V?5y  +  vnny 


(96) 


and 
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For  u  and  v  equal  to  zero  on  the  surface,  Eqs  (96)  and  (97) 
can  be  combined  to  yield 
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Now,  since  rj  =  9  corresponds  to  the  airfoil  surface,  this 
rj-line  is  a  streamline  along  which  ip  =  S  everywhere. 
Consequently,  the  rate  of  change  of  ip  with  respect  to  y  , 
must  be  zero  along  this  line.  Then,  noting  that  the  metrics 


are  not  zero  in  general,  the  requirement  for  zero  velocity 
on  the  airfoil  surface  can  only  be  met  by  prescribing 


e^=8  at  q  =  8  (99) 

which  is  the  implicit  form  of  the  surface  boundary  condition 
for  y.  The  vorticity  on  the  airfoil  surface  can  be 
determined  from  the  definition  of  vorticity, 
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Since  v  =  8 
of  change  of 
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everywhere  on  the  surface,  however,  the  rate 
v  with  respect  to  x  is  also  zero.  So,  the 
the  surface  simplifies  to 
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This  boundary  condition  can  be  written  in  terms  of  the 
computational  domain  variables  by  applying  the  chain  rule, 
which  gives 


"  =  5yH5?(Vv  *  Vr)  *  ny-!j(Vy  *  Vy) 


(  1®2) 


This  equation  can  be  simplified,  however,  by  noting  that 
both  and  the  rate  of  change  of  tp^  with  respect  to  |  are 
zero  at  r?  =  8  ,  so  that 


"  w  =  vS?(Vy)  at  n  =  8  {1fi3) 

which  is  the  implicit  form  of  the  surface  boundary  condition 
for  u. 

The  boundary  conditions  in  the  far  field  are  derived  by 

assuming  the  outer  boundary  is  far  enough  away  from  the 

airfoil  that  uniform  flow  at  free  stream  conditions 

prevails.  Thus,  for  £  =  8  and  5  =  |  , 
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The  boundary  conditions  along  the  wake  cut  depend  on 
the  airfoil  and  angle  of  attack.  For  a  symmetric  airfoil  at 
zero  angle  of  attack,  the  wake  cut  can  be  treated  as  a 
streamline  and  the  appropriate  boundary  conditions  are 


tp  =  8 
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(  186  ) 


These  boundary  conditions  assume  uniform  flow  in  the 
streamwise  direction  and  no  component  of  velocity  in  the  y 


or  17  direction.  Therefore,  the  wake  boundary  condition  must 
be  modified  for  non-symmetr ic  airfoils  or  airfoils  at  angle 
of  attack.  Hafeez  (27)  is  currently  working  on  this 
mod i f i cat i on . 

Mod i f i cat i on  0 f  the  Original  Code .  The  original 
computer  program  solved  for  flow  over  a  symmetric  airfoil  at 
zero  angle  of  attack.  The  code  was  written  to  take 
advantage  of  the  symmetry  by  computing  the  solution  for  only 
half  the  domain.  As  a  first  step  in  generalizing  the 
program  to  non-symmet r i c  airfoils,  the  author  and  fellow 
student  Capt  Faran  Hafeez,  modified  the  original  code  to 
extend  the  computations  over  the  entire  domain.  The 
modified  code  was  tested  by  applying  it  to  a  NACA  Q8  1  2 
airfoil  at  zero  angle  of  attack,  for  which  the  solution  had 
been  previously  computed  with  the  original  code.  Due  to 
problems  in  applying  the  implicit  boundary  conditions 
correctly,  all  the  time  available  to  devote  to  this 
algorithm  was  spent  extending  the  scheme  to  the  full  domain 
and  testing  the  modifications.  Therefore,  two  further 
modifications  of  this  code  are  needed  before  it  can  be 
applied  to  an  iced  airfoil. 

First  of  all,  a  suitable  wake  boundary  condition  must 
be  derived  for  flows  at  angle  of  attack  or  flow  over 
non-symme t r i c  airfoils.  One  possibility  is  to  average  the 
velocities  across  the  wake  cut,  as  is  done  in  the 
Beam-Warming  code.  The  task  is  to  obtain  expressions  for 
the  average  velocities  in  terms  of  u  and  however.  As 
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mentioned  previously,  Hafeez  (27)  is  working  to  incorporate 
a  suitable  wake  boundary  condition.  The  second  modification 
needed  is  the  addition  of  a  turbulence  model.  The  separated 
flow  behind  the  glaze  ice  shape  on  an  iced  airfoil  is 
definitely  turbulent  and  accurate  results  cannot  be  expected 
unless  this  phenomena  is  modeled.  Hafeez  (27)  is  in  the 
process  of  modifying  the  code  to  handle  turbulent  flow. 
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IV.  Results  and  Discussion 


The  Beam-Warming  algorithm  was  applied  to  an  iced  NACA 
•912  airfoil  at  angles  of  attack  ranging  from  2-8  degrees. 
The  criteria  used  to  define  a  fully  converged  or 
steady-state  solution  are  presented  in  this  section.  Also, 
the  global  results  of  this  study  are  presented  in  the  form 
of  lift  and  drag  curves,  and  pressure  coefficient  plots  for 
the  flowfield.  Some  local  flowfield  results  are  presented 
as  well.  These  include  u-velocity  component  contours  over 
the  leading  edge  region  for  2.8,  4.8,  6.8,  and  8.8  degrees 
angle  of  attack,  and  velocity  vector  plots  covering  the 
separated  flow  regions.  The  u-velocity  contours  illustrate 
the  separation  and  reattachment  of  the  flow  over  the  ice 
horns,  and  the  velocity  vector  plots  detail  the  reverse  flow 
in  the  separation  bubbles.  The  velocity  profile  on  the 
upper  surface  at  x/c=8.8  is  also  plotted  and  compared  to 
Spring’s  experimental  data.  Computed  and  experimental  flow 
reattachment  locations  are  compared  as  well. 

In  a  separate  section,  the  results  of  the  work  done  to 
extend  Halim’s  formulation  to  the  full  domain  are  presented. 
Streamline  contours  for  a  clean  NACA  8812  airfoil  at  zero 
angle  of  attack  are  presented,  along  with  a  plot  of  the  skin 


friction  coefficient  over  the  airfoil. 
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The  Beam-Warming  code  developed  by  Visbal  was  used  to 
compute  the  fiowfield  over  an  iced  NACA  8812  airfoil  with 
the  same  geometry  as  that  tested  experimentally  by  Bragg  and 
Coirier.  Computations  are  for  flow  at  a  free  stream  Mach 

C 

number  of  8.12  with  a  Reynolds  number  of  1.41X18  . 

Transition  from  laminar  to  turbulent  flow  is  fixed  at 
x/c=~8 .818  on  the  upper  surface,  and  at  x/c=.8327  on  the 
lower  surface.  These  locations  correspond  to  grid  points 
one  step  upstream  of  flow  separation.  Steady-state 
numerical  solutions  were  obtained  for  angles  of  attack  from 
2-8  degrees.  Comparison  of  these  solutions  is  made  with 
experimental  data  and  other  numerical  results. 

The  grid  used  in  these  computations  is  comprised  of  299 
points  in  the  direction  around  the  airfoil,  and  65  points 
normal  to  the  airfoil.  There  are  281  points  on  the  airfoil 
surface  to  define  the  airfoil  and  ice  shape,  leaving  49 
points  on  each  side  of  the  wake  cut  extending  to  the  outer 
boundary.  The  outer  boundary  is  located  28  chord  lengths 
from  the  airfoil.  The  minimum  spacing  normal  to  the  airfoil 
surface  is  2.8X18  chord  lengths.  The  minimum  spacing 
along  the  airfoil  surface  is  .881  chord  lengths,  specified 
at  the  airfoil  trailing  edge.  Deviation  from  the  normal 
NACA  8812  shape  to  describe  the  ice  geometry  begins  on  the 
upper  and  lower  surfaces  at  14%  chord.  At  these  positions, 
a  minimum  spacing  along  the  surface  of  .82  chord  lengths  was 
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specified.  The  ice  shape  geometric  data  used  are  those 
reported  by  Khodadoust  (11).  A  close-up  of  the  leading  edge 
portion  of  his  grid  is  shown  in  Figure  3. 

Convergence  Criteria.  Before  discussing  the 

steady-state  results  obtained  with  the  Beam-Warming  code,  it 

is  necessary  to  specify  how  the  steady-state  condition  was 

defined.  Visbal  (23)  has  shown  that  convergence  can  be 

determined  by  monitoring  the  lift  and  drag  coefficients  as 

the  solution  progresses.  This  procedure  was  used  to  define 

convergence  in  the  current  study.  A  solution  was  considered 

fully  converged  when  the  lift  coefficient  changed  by  less 

than  .1%,  and  the  change  in  drag  coefficient  was  less  than 

.  8882  over  180(9  iterations.  When  these  conditions  were  met, 

-  8 

the  root  mean  square  of  the  residual  was  less  than  1.8X18 
for  every  solution  obtained. 

To  illustrate  the  convergence  behavior  for  a  typical 
case,  the  root  mean  square  of  the  residual,  and  the  lift  and 
drag  coefficients  for  a=4.8  degrees  are  plotted  versus  the 
number  of  iterations  in  Figures  4-6,  respectively.  This 
solution  was  obtained  using  the  solution  for  oc=2.8  degrees 
as  an  initial  condition.  With  a  previous  solution  as  an 
initial  condition,  note  that  the  oscillations  in  the  lift 
and  drag  coefficients  damp  out  very  quickly.  Similar 
results  were  obtained  for  a=6.8  and  a =8.0  degrees. 
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Globa  I  Resu 1 t  s .  The  computed  values  of  airfoil  lift 
coefficient  versus  angle  of  attack  are  presented  in  Figure 
7.  Also  plotted  on  the  figure  are  the  experimental  data  of 
Bragg  and  Coirier,  and  the  numerical  results  of  Potapczuk 
and  Cebeci.  The  lift  curve  for  a  clean  airfoil  is  included 
as  well  to  illustrate  the  overall  effect  of  ice  on  the  lift. 

The  computed  results  from  the  Beam-Warming  code  compare 
reasonably  well  with  the  experimental  data  for  angles  of 
attack  from  2-6  degrees,  with  a  maximum  difference  of  11.5% 
occurring  at  a=6  degrees.  There  is  an  obvious  discrepancy 
between  the  experimental  data  and  the  computed  solution  for 
a=8.0  degrees,  however. 

The  experimental  data  indicate  that  the  airfoil  stalls 
at  an  angle  of  attack  of  7.9  degrees.  Increasing  the  angle 
of  attack  to  8.Q  degrees  results  in  a  4.6%  decrease  in  the 
lift  measured  exper imenta I  I y .  This  decrease  is  not  seen  in 
the  present  numerical  results,  as  the  computed  lift 
coefficient  for  <*=8.9  degrees  is  7.6%  higher  than  the 
experimental  value.  The  reason  for  this  discrepancy  is  that 
the  flowfield  characteristic  is  completely  different  after 
the  airfoil  stalls. 

For  angles  of  attack  below  stall,  the  flow  separates 
over  the  ice  shape  and  then  reattaches  at  some  downstream 
location  on  the  upper  and  lower  surfaces.  Thus,  regions  of 
separated  flow,  or  separation  bubbles,  form  on  the  airfoil. 
At  angles  of  attack  above  stall,  the  flow  never  reattaches 
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to  the  upper  surface,  however.  This  phenomena  is  pointed 
out  in  Khodadoust’s  flow  visualization  study,  and  Potapczuk 
has  noted  this  trend  in  his  numerical  studies  as  well.  This 
same  flow  behavior  is  seen  in  the  present  study  also,  as 
will  be  shown  later 

Therefore,  the  flow  over  the  iced  airfoil  behaves  quite 
differently  before  stall  than  it  does  after  stall.  The 
results  in  Figure  7  indicate  that  the  Beam-Warming  code  does 
a  good  job  modeling  the  separ a ted/r ea t tached  flow  which 
occurs  before  stall,  but  tends  to  over  predict  the  lift  in 
the  completely  separated  flowfield  which  results  after 
sta  I  I  . 

Further  examination  of  Figure  7  shows  that  Potapczuk’s 
result  for  «=8.Q  degrees  correctly  predicts  the  decreased 
lift  expected  after  stall.  To  compute  the  correct 
post-stall  lift,  Potapczuk  investigated  the  unsteady 
behavior  of  the  flow  for  a=8.G  degrees,  which  he 
characterizes  as  periodic  vortex  shedding.  He  obtained  the 
result  plotted  in  Figure  7  for  this  angle  of  attack  by 
averaging  the  pressure  coefficient  values  at  each  station 
over  several  vortex  shedding  periods  (15).  Therefore,  even 
though  a  steady-state  solution  can  be  obtained  for  a=8.Q 
degrees,  it  appears  to  be  an  inaccurate  representation  of 
the  actual  flowfield,  and  results  in  a  computed  lift 
coefficient  that  is  too  high. 


The  drag  curve  for  the  iced  airfoil  is  shown  in  Figure 
8.  Drag  values  computed  in  the  present  study  compare 
favorably  with  the  experimental  data  for  angles  of  attack 
from  2-6  degrees.  The  computed  Cd  differs  from  the 
experimental  value  by  24%  for  a=2.8  degrees,  improving  to 
less  than  18%  difference  at  a=6.8  degrees.  Similar 
results  can  be  seen  for  the  other  numerical  solutions.  AM 
three  schemes  over  estimate  the  drag  slightly  for  small 
angles  of  attack,  while  under  estimating  it  for  larger 
angles  of  attack.  For  a=8.G  degrees,  once  again  the  large 
discrepancy  between  the  numerical  result  from  the 
Beam-Warming  code  and  the  experimental  data  is  attributable 
to  the  unsteady  nature  of  the  flowfield  at  this  angle  of 
attack.  The  steady-state  solution  computed  in  this  study  is 
just  not  accurately  describing  the  post-stall  flow 
character i st i cs . 

The  computed  pressure  coefficient  over  the  iced  airfoil 
is  shown  in  Figures  9-11  for  angles  of  attack  of  2.8,  4.8, 
and  6.8  degrees,  respectively.  Spring’s  experimental  data 
are  plotted  as  well  for  comparison.  In  addition  to  the 
experimental  data,  comparison  is  also  made  with  other 
numerical  results.  For  a=2.8  degrees,  Potapczuk’s 
numerical  results  are  shown,  while  the  numerical  results  of 
Cebeci  are  included  in  the  plots  for  «=4.8  and  a=6.8 
degrees.  It  should  be  noted  that  the  experimental  data  and 
previous  numerical  results  were  obtained  by  scaling  plots 
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contained  in  the  references  previously  cited  for  these 
authors.  Care  was  taken  to  obtain  as  much  accuracy  as 
possible,  but  there  is  certainly  some  error  involved  in 
obtaining  data  by  scaling.  Therefore,  the  intent  of  the 
comparisons  shown  in  these  plots  is  simply  to  reveal  general 
trends  between  previous  work  and  the  current  study. 

Note  that  the  computed  results  of  the  present  study 
show  sharp  peaks  in  the  pressure  coefficient  as  the  flow 
goes  over  the  ice  shape.  These  peaks  are  due  to  the  very 
large  gradients  in  the  velocity  and  pressure  in  this  region. 
The  oscillations  in  the  Cp  profiles  near  the  ice  shape 
result  as  the  code  attempts  to  resolve  the  gradients.  The 
original  data  of  Cebeci  and  Potapczuk  also  showed  sharp 
peaks  in  the  pressure  coefficient  near  the  ice  shape,  but 
these  are  not  seen  in  Figures  9-11  because  too  few  points 
are  shown.  The  experimental  data  do  not  exhibit  the  sharp 
peaks,  however.  This  is  because  the  pressure  taps  do  not 
sense  the  i ns  tan taneous  pressure,  but  rather  some  average 
value  in  the  region. 

Examining  Figures  9-11  shows  that  the  computed  pressure 
coefficient  on  the  lower  surface  agrees  well  with  the 
experimental  data  for  all  angles  of  attack.  Similar 
agreement  can  be  seen  for  the  other  numerical  results  as 
well.  On  the  upper  surface,  however,  agreement  with  the 
experimental  data  is  a  function  of  angle  of  attack.  For 
a=2.Q  degrees,  the  computed  pressures  agree  well  except  in 


the  region  x/c=.18-.28  At  this  location,  there  is  a 
change  in  the  slope  of  the  experimental  data  that  the  code 
does  not  pick  up.  The  numerical  results  are  relatively  flat 
in  this  region.  Potapczuk's  results  agree  more  closely  with 
the  experimental  data  in  the  region  just  behind  the  ice  horn 
on  the  upper  surface. 

The  discrepancy  behind  the  ice  shape  on  the  upper 
surface  is  worse  as  the  angle  of  attack  increases,  as  can  be 
seen  in  Figures  18  and  11.  For  «=4.8  degrees,  the 
experimental  and  computed  results  differ  substantially  over 
about  the  first  28%  chord.  In  this  region,  the  computed 
pressure  coefficients  are  much  higher  than  the  experimental 
results.  The  high  velocity,  low  pressure  peak  in  the 
experimental  data  at  approx ima te I y  x/cs.15  is  not 
predicted  in  the  numerical  results  at  all.  Note  that 
Cebeci’s  data  does  not  compare  well  with  the  experimental 
data  in  this  region,  either.  For  a=6.8  degrees,  the  same 
trend  can  be  seen.  Again,  Cebeci’s  data  and  the  results  of 
the  present  study  are  essentially  the  same,  with  both  codes 
predicitng  pressure  coefficients  much  higher  than  the 
experimental  data.  However,  at  this  angle  of  attack  the 
discrepancy  in  the  results  extends  to  about  48%  chord. 

Loca I  Resu Its.  Local  f I ow field  results  are  included  to 
show  the  size  and  location  of  the  separation  bubbles. 
Separation  over  the  ice  horns  is  illustrated  in  the 
u-velocity  contours  plotted  in  Figures  12-19,  corresponding 
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to  angles  of  attack  of  2.8,  4.8,  6.8,  and  8.8  degrees.  For 
each  angle  of  attack,  one  labeled  and  one  unlabeled  contour 
plot  is  shown.  The  growth  of  the  upper  surface  separation 
bubble  with  increasing  angle  of  attack  can  be  seen  clearly 
in  these  figures.  For  a=6.8  and  a=8.8  degrees,  it  can 
be  seen  from  Figures  16-19  that  the  flow  does  not  reattach 
like  it  did  for  lesser  angles  of  attack.  The  u-velocity 
contours  over  the  entire  airfoil  for  a=6.8  degrees  are 
shown  in  Figure  28.  This  figure  indicates  a  region  of 
reverse  flow  near  the  surface  along  the  entire  airfoil  upper 
surface.  The  data  for  this  angle  of  attack  indicate  the 
flow  does  reattach  at  x/c=.9899.  Figure  21  shows  than  an 
even  thicker  separated  flow  region  extends  over  the  airfoil 
at  an  angle  of  attack  of  8.8  degrees.  At  this  angle  of 
attack,  the  flow  never  reattaches  to  the  airfoil  upper 
sur  face . 

Another  representation  of  the  separation  bubbles  is 
given  in  the  vector  velocity  plots  of  Figures  22-27.  These 
figures  clearly  indicate  the  separated  flow  regions  caused 
by  the  ice  shape,  as  well  as  the  reverse  flow  within  the 
bubbles.  Again,  for  a=6.Q  and  a=8.Q  degrees,  plots  of 
the  entire  airfoil  are  included  to  show  the  separated  flow 
over  the  upper  surface. 

The  separation  and  reattachment  points  for  the 
numerical  results  are  shown  in  Table  1  for  each  angle  of 
attack.  In  all  cases,  the  flow  separates  at  the  same 
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location  on  the  ice  horns,  while  reattachment  location  is  a 
f unc tion  of  angle  of  attack.  As  indicated  previously,  the 
flow  over  the  upper  surface  does  not  reattach  at  all  for 
a=8.G  degrees.  Table  1  shows  that  the  flow  is  essentially 
separated  over  the  entire  airfoil  upper  surface  for  6.S 
degrees  angle  of  attack  as  well. 

Comparison  of  reattachment  locations  with  the 
experimental  data  for  a=2.8  and  a=4.8  degrees  is  shown 
in  Figure  28.  It  can  be  seen  from  this  figure  that  the 
computed  results  exhibit  the  correct  trend  for  flow 
reattachment.  On  the  upper  surface,  the  reattachment  point 
shifts  rearward  as  the  angle  of  attack  is  increased,  while 
on  the  lower  surface  reattachment  shifts  forward  as  angle  of 
attack  increases.  Beyond  this  general  trend,  however, 
comparison  with  the  experimental  data  is  poor. 

The  code  predicts  subs t an t i a  I  I y  smaller  separation 
bubbles  on  both  the  upper  and  lower  surfaces.  On  the  upper 
surface,  the  computed  reattachment  location  is  24%  less  than 
the  experimental  value  for  o=2.8  degrees,  and  31%  less  at 
a=4.G  degrees.  On  the  lower  surface,  the  computed  result 
is  approximately  48%  less  than  the  experimental  value  for 
both  angles  of  attack.  If  the  experimental  uncertainty  is 
considered,  the  comparison  is  a  little  better  but  the 
differences  are  still  substantial. 

The  computed  separation  bubbles  are  not  only  shorter 
than  the  experimental  ones,  but  they  are  much  thinner  also. 
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Figure  28.  Comparison  of  Computed  Reattachment  Locations 

with  Experiment  (Experimental  Values  from  18:52) 


MB25SZZ i 


V-  - 


This  can  be  seen  from  the  velocity  profiles  plotted  in 
Figure  29.  This  plot  compares  the  computed  velocity  profile 
on  the  upper  surface  at  x/c=9.9  to  Spring’s  experimental 
data  for  an  angle  of  attack  of  4.9  degrees.  The  two 
velocity  profiles  are  considerably  different.  The 
experimental  data  show  reverse  flow  in  the  bubble  from  the 
surface  up  to  approximately  y=.5  inches.  Computed  results 
indicate  reverse  flow  only  up  to  approximately  y  = . 1 5 
inches.  In  general,  the  computed  bubble  is  roughly  half  the 
thickness  of  that  measured  experimentally. 

This  mismatch  between  computed  and  experimental 
velocity  profiles  at  a  location  near  the  ice  shape  is  not 
unique  to  the  present  study.  Cebec i  has  also  noted 
discrepancies  between  the  experimental  velocity  profiles  and 
those  computed  with  his  interactive  boundary  layer  technique 
for  locations  within  the  separation  bubble  (29). 


Resu I ts  for  a  C I ean  NACA  99 1 2  Us i ng  Hal i m ’ s  Formu I  at i on 


The  code  developed  by  Halim  was  extended  to  the  full 
domain  during  the  course  of  this  study,  and  the  original 
intention  was  to  perform  the  other  modifications  required  to 
apply  this  code  to  an  iced  airfoil.  The  modifications 
needed  could  not  be  completed  within  the  time  allowed, 
however,  so  the  code  was  never  applied  to  an  iced  airfoil. 
Nonetheless,  a  test  case  for  a  clean  NACA  9912  was  run  with 
this  code  to  check  out  the  modifications  which  were  made  to 
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extend  the  code  to  the  full  domain.  For  completeness,  the 
results  of  this  test  case  are  presented  here. 


The  code  was  run  for  a  NACA  8812  airfoil  at  zero  angle 
of  attack,  with  a  Reynolds  number  of  12,588.  The  salient 
feature  of  this  flow  is  the  separation  that  occurs  near  the 
trailing  edge.  In  Halim’s  previous  study  with  this  airfoil, 
separation  was  predicted  at  x/c=.8178  .  For  the  grid  used 

in  the  present  study,  the  predicted  separation  point  was 
x/c=.9817  Streamline  contours  over  the  airfoil  showing 

the  separated  flow  regions  can  be  seen  in  Figure  38.  in 
Figure  31,  the  skin  friction  is  plotted  versus  x/c.  Halim’s 
(19)  results  for  this  airfoil  are  shown  for  comparison. 

From  this  figure  it  can  be  seen  that  the  skin  friction  goes 
to  zero  near  x/c=8.9. 
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V.  Conclusions  and  Recommendations 


in  the  introduction  to  this  thesis,  it  was  pointed  out 
that  one  of  the  primary  reasons  for  applying  computational 
fluid  dynamics  codes  to  the  iced  airfoil  problem  is  to 
develop  and  validate  a  code  that  can  be  used  as  a  design 
tool.  Such  a  code  could  then  be  used  to  predict  the 
sensitivity  of  proposed  airfoil  designs  to  surface  roughness 
and  icing.  It  is  evident  from  the  results  presented  here, 
however,  that  much  more  research  will  be  needed  before  an 
accurate  computational  design  tool  can  be  fully  developed. 
Clearly,  this  goal  has  not  been  achieved  in  the  present 
study . 

What  has  been  shown  in  this  work  is  that  the 
Beam-Warming  code,  as  implemented,  produces  very  good 
estimates  of  the  airfoil  lift  and  drag  for  angles  of  attack 
below  stall.  For  flow  at  a  free  stream  Mach  number  of  8.12 

A 

and  a  Reynolds  number  of  1.41X18  ,  the  computed  results 
compared  favorably  to  experimental  data  and  previous 
numerical  results.  The  situation  is  quite  different  for 
angles  of  attack  above  stall,  however.  Although  a  converged 
steady-state  solution  was  obtained  for  a=8.8  degrees,  the 
lift  and  drag  coefficients  obtained  from  this  solution  did 
not  compare  well  with  experimental  data.  The  predicted  lift 
for  this  condition  was  too  high,  while  the  drag  coefficient 


was  too  low  compared  to  experimental  values.  Potapczuk  has 
shown  that  the  flow  is  really  unsteady  at  this  angle  of 
attack,  and  pressures  must  be  averaged  over  several  vortex 
shedding  periods  to  obtain  an  accurate  estimate  of  the  lift 
and  drag. 

While  some  success  was  achieved  in  the  computation  of 
the  global  airfoil  performance  parameters,  the  computed 
local  flowfield  results  generally  compared  poorly  with  the 
experimental  data.  In  particular,  the  separated  flow  behind 
the  ice  horns  reattached  sooner  than  expected  for  «=2.8 
and  a=4.Q  degrees,  resulting  in  computed  separation 
bubbles  much  smaller  than  those  measured  experimentally.  A 
comparison  of  the  experimental  and  computed  velocity 
profiles  in  the  upper  surface  separation  bubble  showed  that 
the  computed  separation  region  is  also  much  thinner  than 
that  measured  experimentally.  Clearly,  the  local  results 
obtained  in  this  study  are  of  questionable  value. 

There  is  no  doubt  that  the  Beam-Warming  implementation 
of  the  full  Nav i er-Stokes  equations  can  be  used  to  obtain 
accurate  flowfield  information  for  the  iced  airfoil  problem, 
however.  Given  that  reasonable  results  have  been  obtained 
by  Cebeci  and  Potapczuk  using  approximate  forms  of  the 
governing  equations,  one  would  certainly  expect  to  achieve 
results  as  good  as  or  better  than  these  from  the  full 


Nav i er -Stokes  equations.  To  do  so,  however,  will  require 
further  investigation  in  at  least  three  areas. 

One  aspect  of  the  problem  requiring  further  study  is 
the  effect  of  the  grid  on  the  accuracy  of  the  results.  In 
this  study,  only  two  grids  were  applied  to  the  problem.  The 
first  of  these  grids  was  used  to  generate  the  solutions 
presented  here,  and  has  been  described  earlier  in  this 
thesis.  A  second  grid  with  about  four  times  as  many  grid 
points  concentrated  around  the  ice  shape  was  also  applied  to 
the  problem.  The  minimum  spacing  normal  to  the  wall  was 
less  by  a  factor  of  IQ,  as  well.  Unfortunately,  a  converged 
solution  could  not  be  obtained  with  this  grid.  This  may 
have  been  due  to  the  fact  that  clustering  more  points  around 
the  leading  edge  left  too  few  points  remaining  over  the  aft 
portion  of  the  airfoil.  Nonetheless,  this  does  illustrate 
the  point  that  the  particular  grid  used  has  a  major  impact 
on  the  resulting  solution.  Given  adequate  time  to  properly 
refine  the  grid  for  the  iced  airfoil,  there  is  every  reason 
to  expect  that  better  results  could  be  obtained. 

The  transition  location  specified  is  another  aspect  of 
the  problem  requiring  further  study.  The  location  of 
transition  has  a  definite  effect  on  the  lift  and  drag 
results.  Computations  were  made  for  «=2.Q  degrees  with 
the  transition  location  fixed  at  x/c=.Q5  on  both  the  upper 
and  lower  surfaces.  Shifting  the  transition  to  thess 


locations,  from  x/c=-.*18  on  the  upper  surface  and 
x/cs.8327  on  the  lower  surface,  resulted  in  a  lift  increase 
of  13.8%,  and  a  decrease  in  the  drag  of  13. 8%.  Based  on 
these  results  for  just  two  different  transition  locations, 
it  certainly  seems  worthwhile  to  systematically  investigate 
the  effects  of  transition  on  both  the  global  and  local 
results  for  a  range  of  transition  points. 

The  final  item  recommended  for  further  study  is  the 
turbulence  model.  Potapczuk  has  pointed  out  that  the 
Ba I dwi n-Lomax  turbulence  model  does  not  always  compute  the 
correct  eddy  viscosity  in  a  massively  separated  turbulent 
boundary  layer.  This  problem  arises  from  the  way  the  length 
scale  in  the  outer  region  of  the  boundary  layer  is 
calculated.  Based  on  the  length  scale,  it  is  possible  for 
the  resulting  eddy  viscosity  to  be  too  high  or  too  low.  in 
the  present  study,  no  systematic  investigation  of  the 
turbulence  model  was  accomplished.  Frankly,  such  an 
undertaking  would  be  a  thesis  topic  in  itself.  The 
recommendation  for  future  study  of  the  model  is  based  solely 
on  Potapczuk’s  reported  experiences  with  the  model.  insofar 
as  he  has  identified  the  model  as  a  potential  problem  in 
computing  the  flow  over  an  iced  airfoil,  investigation  of 
this  model  would  definitely  be  an  important  part  of  any 
follow-on  effort  designed  to  improve  upon  the  results  of  the 


present  study. 


Appendix  A:  Nondimensional i 2 a t i on  0 f  the  Govern i ng 
Equat i ons 

The  procedure  used  here  to  normalize  the  governing 
equations  is  based  on  the  method  outlined  in  Anderson  and 
Tannehill  (21).  As  shown  in  that  refernece,  all  terms  in 
the  Nav i er -Stokes  equations  can  be  put  in  dimensionless  form 
by  substituting  the  following  definitions  into  the 
dimensional  form  of  the  equations: 

*  *  *  t 

X  =  x/c  y  =  y/c  t  = 

CD 

*  ★  * 

u  =  u/V  v  =  v/V  p  =  p/p 

CD  CD  r*  l* 

u  -  e*  =  e/M*  T*  =  T/T^  <  A 1  ) 


The  procedure  is  illustrated  for  representat i ve  terms.  One 
shear  stress  term  is  non-d imens i ona I i zed ,  along  with  one 
heat  flux  term,  the  total  energy  term  and  the  equation  of 
state. 

Shear  St  ress  Term 

In  dimensional  form  the  shear  stress  term,  t  , i s 

xx 

g i ven  by 


Txx  =  21  «  *  £  lux  -  2/3  (  ux  *  vy  ) 


( A2  ) 
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Substituting  from  Eq  (A1)  for  each  term  in  Eq  (A2)  yields 


afu*V  ] 

T*« =  H  v  ♦  v  J  -r— 

T  , 


r  *  *  'i 

r  a 

r  * 

u  V 

l  « 

3 

_  +  _ 

r  * 

V  V 

00 

V 

r 

{  W  *  H*e  J 

.  9 

a 

(-) 

- 

Factoring  out  the  constants  u  ,  V  ,  and  c  gives 

CO  CD  ” 


fj  V 
*00  00 

:xx  "  ~ 


{  2(  * £*) 
-  2/3  fj* 


f  it  ^ 

au  av 

2/3  |i  +  e  - 7  +  - * 

^  >  ax  ay 


The  Reynolds  number  based  on  the  chord,  c,  is  defined  as 


p  V  c 
00  00 


which  can  be  rewritten  as 


M  p  V 

'w  0D 


so  that 


( A7  ) 


u  V 
K®  ® 


P  v2 
K®  ® 

-RTi- 


Comparing  Eq  (A7)  to  Eq  (A4),  it  can  be  seen  that  txx  can  be 

2 

made  dimensionless  by  dividing  by  pv^.  Therefore, 


xx 


Pv 


xx 

T 


® 


(A8) 


or 


★ 

+  e  ) 


★  * 

-  2/3  (p  +  e 


av 

i 

ay 


<A9) 


The  dimensionless  expressions 
in  the  same  manner . 


for 


■xy 


and  z 


yy 


can  be  der i ved 


Heat  F I ux  Term 

The  dimensional  heat  flux  term  in  the  x-direction  is 


d 


x 


(  A  1  8  ) 


Substituting  for  the  dimensional  variables  gives 


92 


which  can  be  rearranged  as 


Now,  c 


and  T 

CO 


Subst  i 
( A 1 2  ) 

d 


q 


x 


C  T  n 

p  OF® 


can  be  written  as 
P 


c 


P 


can  be  expressed  in  terms  of  the  Mach  number  as 


T 

® 


OD 


tuting  Eqs  (A5),  (A6),  (A13),  and  (A14)  into  Eq 
yields 


=  P  V3 

00 


-  1 

r  * 

P  + 

*  1 

e 

3T*  ) 

( y-  1  )  M2  Re 
•  ®  c 

P'T  j 

ax*  j 

from  which 


Simi lar ly, 


Total  Energy  Term 

The  total  energy  is  given  by 


E 


t 


=  P 


r 

e  ♦ 


u2  ♦ 


J 


Substituting  for  the  dimensional  variables  results  in 


E 


t 


P  v2 

CD  CD 


I  *  *  * 

p  e  +  p 


from  which 


of  State 


Equat i on 

Dimensionless  forms  of  the  pressure  and  temperature  can  be 
derived  from  the  equation  of  state  given  by 

p  =  pRT  (A22) 

The  gas  constant,  fl,  is  related  to  the  specific  heat  at 
constant  volume,  cy,  by 

R  =  c v { y- 1 )  ( A23 ) 

and  cy  is  in  turn  related  to  the  temperature  by 

e  =  cyT  (A24) 

Combining  Eqs  (A23)  and  (A24)  with  Eq  (A22)  gives 

P  =  pe{y-1)  (A25) 

Substituting  for  p  and  e  from  Eq  (A1)  yields 

P  =  p\v*  e*(y-1)  (A26) 


which  shows 


* 

P 


P 


P 


00  00 


(A27) 


or 
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*  *  * 

P  =  p  e  (7-1 ) 


The  temperature  can  be  non-d imens i ona I i zed  from 


by  substituting  for  the  dimensional  variables  as 


*  2 

*  P  P  v 

V  =  — r— 

p  p  R 


But  V  is  related  to  Mach  number  by 
® 


2  2 

V  =  M  vRT 
®  ®  *  ® 


Making  this  substitution  into  Eq  ( A3G )  gives  the 


dimensionless  temperature 


2  * 
*  *M®  P 

|  S  '  J.  . . 


With  these  definitions,  all  variables  have  been  put  in 


dimensionless  form. 


Appendix  B:  Trans  forma  t i on  to  Genera  I i zed  Coordinates 

The  transformation  of  the  governing  equations  to 
generalized  coordinates  presented  in  this  appendix  is  based 
on  the  method  outlined  in  Anderson  and  Tannehill  (21).  Prom 
that  reference,  the  general  transf ormat i on  from  the  physical 
domain  (x,y)  to  the  computational  domain  (£,q)  is  given  by 


(B1 ) 
( B2  ) 

This  transformation  is  accomplished  by  applying  the 
following  chain  rule  to  the  Navier  Stokes  equations: 


5  =  €(*,y) 
n  =  n(x.y) 


♦ 


n 


an 

x  dr) 


( B3 ) 


3(  •  ) 
ay 


a(-) 


•y  as 


ar 

drj 


(84) 


where  £  ,  n  ,  £  ,  and  n  are  the  metrics  of  the 
x  x  y  y 

transformation.  Differential  changes  in  the  x  and  y 
directions  in  the  physical  domain  are  related  to 
corresponding  changes  in  the  computational  domain  by 


d5  =  €xdx  +  £ydy  ( B5 ) 

drj  =  r)xdx  +  r)ydy  (B6) 


97 


Rewritten  in  matrix  form,  these  equations  are 


( 


A  similar  expression  holds  for  the  inverse  transformation 


from  (f-,n)  to 

( x , y ) ,  which  can 

be  written  as 

dx  = 

*  V" 

( 

dy  = 

y{d5 

*  y 

( 

or  in  matrix 

form 

dx 

dy 

*{ 

.  y5 

X 

n 

yn  J 

dr? 

(B 

Substituting  for  [dx,dy]  in  Eq  (B7)  from  Eq  (BtQ)  gives 


?x  *y 

1 - 

cr 

X 

u/ 

X 

d£ 

drj 

"x  V 

.  y5  yn  . 

dt) 

which  means  the  transformation  metrics  are  related  by 


Recall  that  for  a  square  matrix  [A] 


( B 1  3 


so  it  foil ows  that 


n 


y 


-n 


X 


X 

n 


Defining  the  determinant  above  as  the  Jacobian  of  the 
transformation,  J,  the  metrics  are  related  by 


5 

5 


x 

y 


where 


=  My  -  My 

With  the  relationship  between  the  transformation 
metrics  defined,  the  chain  rule  can  be  applied  to  the 
Nav i er-S tokes  equations  given  by 


3U  3E  x  9F  . 

3T  +  a?  +  5y  - 


8 


Applying  the  chain  rule  yields 


It  is  cosirable  to  have  this  equation  in  strong  conservation 
ft  form.  Anderson  and  Tannehill  (21:254)  suggest  this  can  be 

done  by  dividing  through  by  the  Jacobian  and  rearranging  into 
conservat i ona I  law  form  by  adding  and  subtracting  like 
ft  terms.  For  example, 


“It  (  “T" )  =  -or  “ir" +  u  “§r  [  )  (822) 


Since  J  is  not  a  function  of  t,  this  reduces  to 


ut  =  J 


m. 


(B23) 


For  the  next  term,  r£.,  dividing  by  J  and  expanding  gives 

*  5 


f«».E  1 

— 

f  E  1 

+  E 

r  f  'I 

l  J  . 

5 

l  J  J 

from  which 


f  E  =  J 


r 


-  E 


r  ' 


(B24) 


( B25  ) 


Using  this  procedure  for  the  remaining  terms  in  Eq  (B21) 
gives  the  following: 


C 


e 


1S1 


which  equals  zero.  Similarly,  the  second  term  in  the  last 
brace  of  Eq  (B29)  reduces  to  zero.  Therefore,  the  strong 
conservative  form  of  the  Navier-Stokes  equations  in 
generalized  coordinates  is 


f  u  1  * 

f  E*x  ♦  F?y  1 

4. 

Er>x  *  F,>y  ' 

l  Jt 

l  J  J 

5 

=  S  ( B34 ) 


or 


U.  ♦  E_  +  F  =8 

t  5  n 


where 


u 

=  U/J 

E  =  1/J 

l  5»E  *  V  ) 

F  =  1/J  j 

.  "x*  *  nyF  ) 

(B35  ) 


( B36  ) 


This  equation  can  be  further  manipulated  to  separate  the 
viscous  and  inviscid  terms,  as  well  as  the  cross  derivative 


terms.  Define  a  vector  E«  to  contain  the  inviscid  flux  and 


pressure  terms  in  E,  and  a  vector  Eg  to  contain  the 

A 

corresponding  terms  in  F.  The  components  of  E1  are 


E 


1 


5XPU  ♦  ^yPV 

Cx(pu2+p)  +  5y  puv 

5XPUV  ♦  5y[pV+P2] 

*x(Et+p)u  +  5y[et+p)v 


which  can  be  rewritten  as 


E 


1  " 


*  5yv) 

PU(5*U  *  «/>)  *  Ex’ 

pv(sxu  -  5yv)  *  V 
(Et*’)  (?xu  ♦  Eyv) 


The  velocity  components  in  the  computational  domain, 
contravar i ant  velocities,  are  defined  as 

U  =  5XU  *  v 

v  5  nxu  ♦  nxv 

Using  these  definitions,  E1  can  be  written  as 


(B37) 


(838) 


termed 

( B39 ) 
(B48) 
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E1  = 


puU  +  £XP 
pvu  +  5yp 

(Et  +  P)  u 


(B41) 


In  a  similar  manner , 


E2  = 


puv  +  nxP 
pvv  +  nyp 

(E.*p)  v 


( B42 ) 


The  remaining  viscous  flux  and  heat  flux  terms  are  separated 
into  four  vectors  such  that  the  cross-derivative  terms  are 
isolated.  These  vectors,  in  general,  are  expressed  as 


V  =  ♦  v2(G,un) 

w  =  .  w2(0.0n) 


(843) 


( B44  ) 


where  and  Vg  represent  the  terms  in  E  and  and  Wg 
represent  terms  in  F.  The  final  form  of  the  transformed 
Nav  i  er -Stokes  equations  are 


'  *  \  +  %  =  +  25l  nJ 


W,  [u , U5)  *  W2  (u , uj  (B45) 


The  components  of  vectors  V2,  W1 ,  and  can  be 
assembled  after  applying  the  chain  rule  to  each  viscous  and 
heat  flux  term--xxx,  xxy,  xyy,  qy,  and  qy.  For  example, 
rxx  is  given  by 


txx  =  2  in+c)  u x  -  2/3  (p+e)  £  u 


+  v 

X  y 


(846 


After  applying  the  chain  rule,  this  equation  becomes 


=  4/3  (p+e)  [5xuc  +  nxun] 

-  2/3  (p+e)  ^yv^  +  Hyvnj 


(B47) 


The  other  viscous  terms  become 


Txy  =  <M^>  (  uE5y  *  unny  *  »,={„  *  vn„x  ) 


( B48  ) 


Tyy  =  4/3  «"*e>  («yve  *  Vn) 


-  2/3  ^xu^  +  nxunj  <B49 


Applying  the  chain  rule  to  the  heat  flux  terms  yields 


qx  =  -  cp  (  -£r  +  "WT  )  (  Vx  +  Vx  )  (B58) 

qy  =  "  cp  (  “ftr  +  ~hT  )  (  Vy  +  Vy  )  (B51) 


Visbal  (23)  has  shown  that  these  terms  can  be  assembled  such 
that  V ^  and  are  functions  of  derivatives  of  g  only, 
while  Vg  and  W2  are  functions  of  derivatives  of  17  only. 

The  result  of  rearranging  the  terms  in  this  manner  follows: 


V1  = 


b  u  ♦  b  v 
1  {  2  g 


b2Ug  *  b3Vg 


b  t  uu^  +  b2  ^  u v ^  +  vu^ 


+  b_vv  +  b  .T 
3  g  4  g 


(B52) 


S 


C 


u  + 

n 


c 


2 


V 

n 


CQU  +  c.v 
3  r)  4  n 


C.UU  +  C-UV  +  c.vu  +  c.vv  + 
1  n  2  n  3  T]  4  rj 


ccT 
5  rj 


(B53) 


C1Ug  +  C3Vg 


C2Ug  +  C4Vg 


( B54 ) 


C.UU„  +  C.VU,  +  C-UV,  +  C.VV,  +  C -T_ 


8 

diun  *  c 

d2un  *  * 

d„uu  +  d„|  uv  +  vu 
1  rj  2L  H  T 

where 

b,  =  (M ♦€>  (  4/3  5*  4  C*  ] 

b2  =  lf«l  [  1/3  5x5y  ) 
b3  =  (M«)  (  ?x  4  4/3  {2  ) 

b4  *  cp  (  “£r  +  )  (  5x  ' 

c1  =  -  (M+e)  ^  4/3  £xnx  +  5yn, 

c g  —  ~  { /j+  e )  ^  2/3  5xny  -  5yn. 

c3  =  (^+e)  (  5xny  '  2/3  5yn 

c4  =  -  (M+O  (  5xnx  +  4/3  5yn 

c5  =  "  (  “f*7  +  “Sr^  )  (  ?xnx 

dt  =  (M+e)  [  4/3  nx  +  ny  j 
d2  =  (n  +  c)  ^  1/3  rjxny  j 
d3  =  ( m+g  >  (  nx  ♦  4/3  ny  ] 

d4  =  cp  (  “£?  +  )  (  ^x 


c 


€ 


Appendix  C:  Jacob i an  Mat r i ces  for  the 
Beam-Warming  Scheme 

The  Jacobian  matrices  which  arise  during  linearization 
of  the  governing  equations  are  derived  in  this  section.  Two 
inviscid  flux  vectors  and  two  viscous  flux  vectors  need  to 
be  linearized.  The  Jacobian  for  one  of  each  will  be  derived 
to  illustrate  the  method,  and  the  other  results  simply 
stated.  Derivation  of  the  linearized  equations  follows  the 
method  outlined  in  Anderson  and  Tannehill  (21),  and  uses 
Visbal’s  (23)  notation  for  the  Jacobian  matrix  elements. 

Jacob i ans  for  the  Inviscid  Terms 

The  inviscid  flux  vectors  are  linearized  as 


EE"  =  [ 

A  j  AUn 

(Cl  ) 

*E2  =  [ 

B  j  AUn 

( C2  ) 

where  the  Jacobians  are 

r  a  l  s  • 

aE^U)" 

( C3 ) 

L  J 

au 

>]  ■■ 

aE2(U)n 

( C4 ) 

au 

The 

Jacobian  for  E1  is  derived  by  expressing 

in  terms  of 

the 

components  of  U, 

where  U 

is  given  by 

€ 
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Mi.  -*».  - 


In  terms  of  the  components  of  U,  the  E1  components  can  be 
written  as 


(« 


E 


+ 


The  pressure  term  appearing  in  Eqs  (C8),  (C9),  and  (CIO) 
must  be  expressed  in  terms  of  U  components  before  E1  can  be 


differentiated.  The 

pressure  is 

g  i  ven 

by 

1 

II 

O. 

[£t ' 

1/2 

P  (  u* 

!  2 
+  V 

)] 

(Cl  1 

which  can  be  rewritten  as 

P  =  (  T  *  1  ) 

1 

-  f  u2 

)] 

( C 1  2 

1u 

1  l  2 

Now  that  all  components  have 

been  expressed 

in  terms 

of 

the 

U  variables,  the  derivative 

of  E 

1  can 

be  taken  wi th 

r  espec  t 

to  U.  The  resu 1 t i ng 

Jacobian,  j 

*]• 

can  be  expressed 

as 

ail 

3  1 2 

al3 

3  1  4 

M  ■ 

a21 

3  22 

3  23 

a24 

( C 1  3 

L  J 

a31 

a32 

a33 

3  34 

a4  1 

3  42 

a43 

3  44 

where 

^  1 
i  i 

for 

i  = 

1,4;  j 

=  1,4 

( C 1  4 

For  example, 


*1,  =  -^r  *  4i7(u2?x  ♦  u3«, )  =  •  <c,5> 

dE 

8 1 2  =  ““  =  |oj(  U25x  +  U3*y  )  =  {C16) 

dE 

*13  =  =  lt^[  U2«*  *  U3«,  }  =  «y  <01,) 

fu;(  Vx  *  U3«y  )  *  S  (C18) 

4 


The  remaining  elements  are  determined  in  the  same  manner. 
The  r esu It  is 


21 

3 

-uU 

+  *5X 

22 

= 

u  - 

(?-2)u5x 

23 

r 

u«y 

-  (7-1)vfx 

24 

3 

(7~ 

^x 

31 

= 

-VU 

♦ 

32 

= 

v*x 

-  (Tf-1)ufy 

33 

— 

U  - 

( Tf“2  )  v£y 

34 

- 

(If- 

1)5y 
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41  = 

(  2*  * 

>Et 

P 

)“ 

42  = 

-] 

l«x  - 

(7-1)uU 

43  = 

-] 

l5y  - 

(y-1 )vU 

44  " 

yU 

where 


U  =  5xu  +  5yv 


0>  =  1/2  (7-I) 


2  2 
i  +  v 


The  Jacobian  for  the  Eg  vector  is  derived  in  the  same 
way.  Eg  is  given  by 


nxP“  ♦  Hypv 

!  nx(pu2  +  p]  +  nypuv 

nxpuv  +  ny[pv+p2j 

nx  fe  ,*p) u  *  ny  [et*p] 


The  components  of  Eg  can  be  expressed  in  terms  of  the  U 


components  as 


where  the  expression  for  the  pressure  is  the  same  as  derived 
previously  in  Eq  (C12). 

Defining  the  Jacobian  for  E 2  as 

b11  b 1 2  b 1 3  b  1 4 
b21  b22  b23  b24 

b31  b  32  b33  b34 

b41  b42  b43  b44 

and  differentiating  Eg  with  respect  to  U  as  before,  results 
in  the  following  components: 


( C29  ) 


b21  =  -uv  *  *1), 


b22  =  v  '  <r-2)onx 

b23  “  UT>y  "  ^-Dvr}x 
b24  =  (’“1,r»x 

b31  =  'vV  +  ♦Hy 
b32  =  VrJ*  '  (^1)ur»y 

b33  =  v  -  (y-2)vny 

b34  = 


41  = 
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Hi 

p 

}’ 
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-  <p 

j 

K  ' 

(y-Duv 

43  = 

1 

-  0 

2 

k  ‘ 

( y- 1 ) vV 

44  = 

y  v 

where  the  Jacobians  are 


( C37  ) 


(C38) 


As  before,  the  Jacobian  matrices  can  be  written  as 
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and 


( C38 ) 


1  1 

S  1  2 
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21 

S22 

S23 

S24 

31 

S32 
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5  34 

41 

S42 

S43 

S44 

( C48  ) 


The  Jacobian  R  J  can  be  obtained  by  getting  an 
expression  for  U^,  rewriting  in  terms  of  the 
components,  and  differentiating  the  resulting  vector  with 


respect  to  U  .  U  is  obtained  by  differentiating  U,  which 
C  c 

resu Its  in 


where  b^  bg,  bg,  and  b^  are  as  defined  in  Appendix  B. 
Rewritten  in  terms  of  the  components,  the  components  of 


V«  can  be  expressed  as 


(C44) 
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1/P 
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6i(vuUi) 
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*  b!(U3-rUl) 

*  ba(ur''ui) 

*  b2[u(U3-vUl]  *  V(U2 
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( C45 ) 
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In  the  V.  component,  T  must  be  cast  in  terms  of  U 
variables  before  can  be  differentiated.  From  the  perfect 
gas  relations, 


e  =  c  T 
v 


Assuming  cy  is  locally  constant, 


is  given  by 


(C48) 


T5  ’ 


e„  /  c 
5  v 


( C49  ) 


where  e^ 
Et  .  For 


must  now  be  expressed  in  terms  of 
Et  given  as 


the  U. 


componen  t  , 


Et  =  P 


e  + 


2  2 
u  ♦  v 


( C5S ) 
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differentiation  with  respect  to  5  yields 


-c  =  4 


e^+uu^+vv^j  +  p^|e+1/2{u2  +  v2)j 


{ C  5  1  ) 


:tP  =  P(' 


e?+uu?+vvj  +  p?(  Et  /  p 


Solving  for  e^  yields 


Et_  P*Ei 


e e  =  — *  -  —T  '  uVvv? 

^  p  P  s 


so  that  can  be  expressed  as 


(C52) 


(C53  ) 
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(C54) 


This  expression  can  be  rewritten  in  terms  of  variables  as 


T5  =  (  U4  -  U1^  -  “(U*  -  UUl)  -  »(°3  -  VUl)  ) 


(C55) 


Differentiation  of  V1  with  respect  to  can  now  be 


performed  to  find  the  components  of 


r«i. 


For  examp  I e , 
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( C  5  7  ) 
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Applying  the  differentiation  to  the  remaining  components  of 
V j  yields  the  following: 


r 

r 

r 

r 


11  = 
12  " 
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components  of  the  Jacobian  S  can  be  expressed  as 
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